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CRITERIA FOR THE OPTIMAL DESIGN 
OF EXPERIMENTAL TESTS 

By George C. Canavos 
Langley Research Center 

SUMMARY 

This report unifies in a concise manner some of the basic concepts that have been 
developed for the problem of finding optimal approximating functions which relate a set 
of controlled variables to a measurable response. The techniques outlined in this report 
have the potential for reducing the amount of testing required in experimental investiga- 
tions. Specifically, two low-order polynomial models are considered as approximations 
to unknown functional relationships. For each model, optimal means of designing experi- 
mental tests are presented which, for a modest number of measurements, yield prediction 
equations that minimize the error of an estimated response anywhere inside a selected 
region of experimentation. Moreover, examples are provided for both models to illus- 
trate their use. Finally, an analysis of a second-order prediction equation is given to 
illustrate ways of determining maximum or minimum responses inside the experimenta- 
tion region. 


INTRODUCTION 

The purpose of this report is to unify some existing results of response -surface 
techniques that are available for the optimal design of experimental tests and to make 
them easily accessible to researchers in various disciplines. In experimental testing, 
it is vital that the test procedure be carefully planned to acquire only those measure- 
ments that are necessary to optimize the end result. Response -surface methods provide 
optimal means of selecting points to observe a dependent variable and to determine an 
approximating equation for an unknown functional form relating the dependent variable 
and a set of controlled variables. Specifically, consider the approximation of the unknown 
relation 

y = f( u l, u 2 ,..., u k) (1) 

with some low-order polynomial inside a desired region defined by the controlled vari- 
ables Up u 2 j . . u^. Moreover, consider an experiment in which the response y is 

to be measured at several selected values of the controlled variables for the purpose of 



making the approximating polynomial form explicit. A basic question which is essentially 
the focal point of this report now arises: With respect to the controlled variables, where 
should the response be measured so that the resulting prediction equation minimizes the 
error between an estimated and an observed response anywhere inside the selected 
region? In other words, an optimal design is defined as one which yields a prediction 
model that minimizes the error of an estimated response. In the discussion of an answer 
to the question, consideration is given to a survey of response -surface methods as applied 
to design -optimization criteria for two polynomial models which have found wide indus- 
trial application (refs. 1 and 2). It is appropriate at this time to note that the term 
"response surface" has risen out of the geometric identification of the problem in that, 
for k s 2 controlled variables, the solution can be characterized by surfaces of con- 
stant response. 

Usually it is mathematically advantageous to restrict the range of each controlled 
variable. Thus, consider the transformation of the natural variables u^, Ug, . . u^ 
to a new set of normalized variables x^, Xg, . . x^. The implementation of this 

variable change is done according to the following scheme: Assume that, for each natural 
variable u^, high and low values are selected such that they serve as boundaries to the 
selected region. Let uj denote the average of these two values and s^ the absolute 
difference between vq and either the high or low value of the ith natural variable. The 
change from uj to x^ is then defined by the relation 

Xi = 1 ■ 1 (i = 1, 2, . . k) (2) 

& i 

As a result, the high and low values of u^ correspond, respectively, to +1 and -1 for xj. 
Moreover, one unit of the normalized variable is equivalent to s^ units of the natural 
variable. Therefore, subsequent discussion is held in terms of the former with the 
inverse change to the latter readily available from equation (2). 

Consideration is now given to a brief discussion of the form of two low -order poly- 
nomials which often serve as good approximations to the unknown function f . One is a 
linear function of the controlled variables whose form is the relation 

k 

y = i3 0 + ^ ^i x i + 6 (3) 

i=l 

where (3^ (i = 0, 1, . . . , k) denotes the unknown coefficients and e is the random 
error in the response y. This first-order model is often useful in studying a given 
response inside a narrow region of the controlled variables where little, if any, curva- 
ture in f is suspected. However, if the desired region is such that the response is not 
likely to be a linear function of the controlled variables, then the approximating function 
should be at least a second-order model given by 
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( 4 ) 


k k k j-1 

7 = ^0 + I ^i x i + I hiH + E I ^j X i X j + 6 

i=l i=l j=2 i=l 


i=l 


k 

The term |3^x^ represents a pure quadratic contribution while the term 
i=l 


k j-1 

^ ^ij x i x j incorporates any cross-product effect. As an example, for k = 2, equa- 
j=2 i=l 

tion (4) reduces to 

y = 0 o + ^l x l + h x 2 + hl x l + hz x 2 + %2 X 1 X 2 + e 


Since a second-order model contains more unknown coefficients than one of first order, 
more measurements of the response are needed to estimate the unknown coefficients. 

Finally, given a set of measurements of the response, it is possible to write either 
one of these models in the matrix form 

y = X£ + e (5) 

where y is a column vector of n responses, X is an nX m (where m = k + 1) 
matrix of normalized -variable values, 0 is a column vector of m unknown param- 
eters, and e_ is the n x 1 vector of random errors. The vector e_ represents the 
change in the response which cannot be accounted for explicitly. This, of course, includes 
changes due to model inaccuracy as well as measurement error. 


SYMBOLS 


B symmetric matrix of quadratic parameters 

k number of controlled variables 

M orthogonal matrix 

m-, fourth mixed moment in matrix x’x 

n number of response measurements 

p^ fourth pure moment in matrix x’x 

s coordinate point of normalized variables 
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s n solution to equation ^ = 0 

u^, Ug, . . ., u k controlled variables 

w vectors of principal axes 

X matrix of normalized variables 

Xp x 2 , . , x k normalized variables 

y vector of measurements 

a rotatability parameter 

0 vector of parameters 

y vector of linear parameters 

e_ vector of random errors 

X 2 , . . A k eigenvalues of B 
0 

a common variance of random errors 

Notation: 

* an estimate 

’ (Prime) transpose of matrix 
(Underlining) vector 

OPTIMAL FIRST -ORDER RESPONSE -SURFACE DESIGNS 

Suppose that the unknown relation given by equation (1) is to be approximated by a 
first-order polynomial in which the jth response is depicted by the expression 

k 

7j = £ 0 + Yj £i X ij + (j = 1, 2, . . n; n > k) (6) 

i=l 
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In matrix form, this set of n linear equations in m unknown coefficients is written as 
y = X0 + £ 


where 



The jth row of X represents the conditions of the controlled variables under which the 
jth measurement of the response is made. 

By using the principle of least squares (refs. 1 and 3), estimates of the components 
of 0 are readily computed by the result 

0 = (X’X) _1 X'y 

where 0 is the vector of estimates. Moreover, the variance -covariance matrix of 0 
is given by the relation 

var .(g) = ct2(x'X) _1 

where a 2 is the unknown common variance of the error vector £ and represents the 
unaccountable change in the response y. Usually a 2 is estimated by the unbiased 
estimate 

„ o y’y - P - xy 

= — zi 

n - m 

The (i,i) element of a 2 (X'X) -1 coincides with the variance of 0 i? while the (i,j) ele- 
ment corresponds to the covariance between 0j and 0j for i # j. Moreover, since 
by definition the variance is a measure of dispersion, then var(fi') is indicative of the 
precision with which each coefficient in equation (6) is being estimated. 

Consider now the prediction of a response at an arbitrary point, inside the desired 
region, whose coordinates in terms of the normalized variables are given by 
( v l’ v 2 , ***’ v k)‘ The estimateci value of y is computed from the relation 

y = v’0 

where 

V* = [1 Vj V 2 ... V k ] 
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By using matrix manipulation, the variance of the predicted value y is given by 

var(y) = v'(x'x) -1 va 2 = v’ var^jv (7) 

Therefore, the conditions that minimize the residual variance of an estimated response 
y are the same as those that maximize the precision of y, where precision is defined 
as the reciprocal of variance. From equation (7), it is clear that the only quantity sub- 
ject to control is the design matrix X whose elements are the values of the controlled 
variables for which the response is to be measured. Hence, the problem of minimizing 
the variance of y is reduced to finding that design matrix X for which the elements 
of var(^) are minimized. For this purpose, the following theorem is given: If the 
design matrix X is selected such that each column of X is mutually orthogonal to 
every other column in X, then the first-order approximation to the response surface is 
optimal in the sense that the residual variance of a predicted response is minimized. In 
other words, the optimal-design matrix X yields a diagonal matrix X’X of the form 


n 0 ... 0 

0 n ... 0 


0 


0 


n 


( 8 ) 


where n is the number of measurements (ref. 4). The proof of this theorem is pro- 
vided in the appendix. 


The inverse of a diagonal matrix such as equation (8) is 
l/n 0 ... 0 

0 l/n ... 0 

(X’X)" 1 = ' 

j • • • 

• • 

l/n 


( 9 ) 


Thus, for a first-order model, the optimal -design configuration results in an expression 
for the error of a predicted response y given by 

|aVn 0 . , , 


var(y) = Jl Vj v 2 ... v^J 


2 /n 


0 a 


0 0 


0 

0 


a2/ n 


( 10 ) 
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Simplification of equation (10) yields the relation 


var(y) = 





( 11 ) 


Equation (11) implies that the optimal first-order design yields an error, for a predicted 


response y, which is a function of r = 
defined region and the point of interest 




the distance between the center of the 


Vi-1 

(vi,v 2 ,...,v k ). 


One class of design configuration that is most useful for fitting first-order models 
and satisfies the optimality criterion is the 2 k factorial design. The factorial arrange- 
ment is accomplished by using as design points all possible combinations of the high and 
low values of the k controlled variables. This, of course, results in a total of n = 2^ 
measurements of the response. Therefore, the 2 k factorial class of orthogonal designs 
is not only optimal in fitting first-order response -surface models but also provides con- 
siderable savings in time and expense as a result of the relatively small number of mea- 
surements needed for a moderate value of k. 


The following hypothetical example is presented to illustrate the procedures for an 
optimal first-order design. Consider the problem of studying the yield of a process as a 
function of agitation speed S, concentration C, and temperature T. It is believed that 
a first-order model such as equation (6) is sufficient to describe the yield inside a region 
bound by agitation speeds of 200 to 400 rpm, concentrations of 2 to 4 percent, and tem- 
peratures of 30° to 50° C. In accordance with the previous discussion, these natural 
variables are transformed to their normalized counterparts by the equations 

_ _ S - 300 
1 100 “ 

x 2 = C - 3 

and 


O 

Since there are k = 3 controlled variables, a total of 2 = 8 measurements of the 
process yield is needed. 
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The design matrix X and the vector of measurements y are given as 


X = 


x 0 

X 1 

x 2 

x 3 


1 

-1 

-1 

-1 


40 

1 

-1 

-1 

1 


54 

1 

-1 

1 

-1 


38 

1 

-1 

1 

1 


56 

1 

1 

-1 

-1 

y = 

52 

1 

1 

-1 

1 


62 

1 

1 

1 

-1 

! 

49 

1 

1 

1 

1 


59 


where the headings Xq, x^, and so forth, indicate the column of values for the corre- 
sponding controlled variable. Since the dot product of any two distinct columns of X 
equals 0, X’X is a diagonal matrix of the form 


X'X = 


8 

0 

0 

0 


0 

8 

0 

0 


0 

0 

8 

0 


0 

0 

0 

8 


Moreover, by matrix multiplication, 

410~ 

34 
-6 
52_ 

As a result, the coefficients of the first-order model are computed to be / 3 q = 51.25, 
j 3 j = 4.25, ~ -0.75, and £3 = 6.50. Therefore, the prediction equation in terms of 

the normalized variables is 

y = 51.25 + 4.25xj - 0 . 75 x 2 + 6.50xg 

This provides an approximating function to use whenever the prediction of yield is needed 
at any arbitrary point interior to the region of experimentation. If, for example, 

S = 250 rpm ^x^ = -|^, C = 3.5 percent ^Xg = and T = 40° C (X 3 = 0 ), the predicted 
yield is y = 48.75 percent. To compute the variance of this estimate y, an estimate of 
a 2 is computed. It turns out to be a 2 = 6.625 and equation (11) is used to compute 
the variance of the predicted yield 



8 



var(y) = (l + 1 + |) = 1.24 


Finally, by using the relations between normalized and natural variables, the approxi- 
mating equation in terms of the latter reduces to 

y = 14.75 + 0.0425S - 0.75C + 0.65T 
OPTIMAL SECOND-ORDER RESPONSE -SURFACE DESIGNS 


Consideration is now given to a discussion of certain useful designs with attractive 
properties for fitting second-order models. It is assumed that a significant degree of 
curvature is suspected in the unknown function y = f( u i> u 2> ,,, > u k) so as to warrant the 
use of a second -order polynomial in which the jth response is represented by 

k k k h-1 

y i = g 0 + I % x ii + I + 2 I ^ihVhi + e i <i * *• 2 n) (l2) 

i=l i=l h=2 i=l 

As in the previous case, equation (12) in matrix notation becomes 


y = X£ + e 

where now 



2 2 

x ll> X 11 x kl x ll x 21> 1 ' •’ x k-l,l x kl 

2 2 

x 12> • * *> x k2 x 12 *• • •> x k2 x 12 x 22’ • * •» x k-l,2 x k2 
2 2 

x l n > • * •> x kn x ln ’ ' * ,,x kn x ln x 2n’ ' * ' ,x k-l,n x kn 


and 



(13) 
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There are exactly (k + 2)(k + l)/2 unknown coefficients in /3 and, thus, the dimensions 
of X are n X (k + 2)(k + l)/2, where n is the total number of measurements. 

Since curvature is suspected, the response y should be measured for at least 
three distinct values of each controlled variable in order to estimate the quadratic effect. 
Based on the usefulness of the 2^ factorial design in the first-order model, it is logical 
to assume that a 3 k factorial design would be appropriate in fitting the quadratic model. 
However, for moderately large values of k, the measurements required for the experi- 
ment become prohibitively numerous. As a result, Box and Wilson (ref. 5) develop what 
they call a class of composite plans which require relatively few measurements and 
achieve certain optimal properties. 

The most frequently used composite plan is the so-called central-composite design 
which is basically a 2^ factorial design augmented by 2k + 1 additional points to allow 
for the estimation of the (k + 2)(k + l)/2 parameters in equation (12). In terms of the 
normalized variables, the additional 2k + 1 points in a central-composite design are as 


follows: 

X 1 x 2 • • • x k 

0 0 ... 0 

-a 0 ... 0 

a 0 ... 0 

0 -a ... 0 

0 a ... 0 

• • • 

• • • 

• # # 

0 0 ... -a 

0 0 ... a 


This design provides five distinct values of each controlled variable for measuring the 
response y if | a | * 1; namely, ±1, ±a, and 0 where (0,0, ...,0) corresponds to the cen- 
ter of the experimentation region. Therefore, in a central-composite design, 2^ mea- 
surements are taken according to a 2^ factorial arrangement, one measurement is 
taken at the center, and the remaining 2k measurements are taken at a distance a 
from the center of the design region. The following sketch provides a geometric repre- 
sentation of the experimentation region for a central-composite plan with three con- 
trolled variables: 
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Depending on the choice of a, the central -composite design is nearly optimal in 
minimizing the error between a predicted and an observed response at any arbitrary 
point inside a defined region for a second-order model. One optimality criterion is a 
property called rotatability. A design is said to be rotatable if the variance of a pre- 
dicted response is unaltered when the design is rotated through an arbitrary angle about 
the center of the experimentation region. To illustrate the effect of rotatability, con- 
sider a central-composite design with two controlled variables and a == 1, The 
design matrix X according to equation (13) is given by 
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and 


(X’X)" 1 = 


5/9 0 

0 1/6 

0 0 

-1/3 0 

-1/3 0 

0 0 


0 -1/3 

0 0 

1/6 0 

0 1/2 

0 0 

0 0 


-1/3 0 

0 0 

0 0 

0 0 

1/2 0 

0 1/4 


Assume the measurements are taken according to this design matrix, and the approxi- 
mating equation is thus determined. 

Consider now the prediction of a response at a point with coordinates (a^ag). The 
estimated value is given by 

y = a’£ 

where 

- “ a l a 2 a l 2 a 2 2 a l a 2^ 

A 

and j3 is the vector of estimates corresponding to equation (12). As before, the vari- 
ance of y is given by the relation 

var (y) = a' (X’X) " 1 aa 2 (1 5) 


After appropriate substitution into and simplification of equation (15), 


var(y) = cr 2 (| - Isij 2 -I 


la , 4 
2 1 


1 a 0 ^ 

2 2 


i a 2 a 2 \ 

4 a l a 2 J 


(16) 


which is a function of location within the region of experimentation. Specifically, con- 
sider the variance of a predicted response at two points which are equidistant from the 
center of the desired region and whose coordinates are given by (0.7, 0.7) and by (\/0.98,o) . 
Substituting into equation (16) and simplifying yield 

var(y)^ ? ^ = 0.3656a 2 


while 

var(y). . — x = 0.5457a 2 

Thus, even though these two points are equidistant from the center, the error in the pre- 
dicted response of the former is obviously smaller than in that of the latter. Hence, a 
design which is not rotatable discriminates against direction. This occurrence is indeed 
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a disadvantage since it is likely that any two points which are equidistant from the experi- 
mental center are equally important. It is apparent that an optimal (orthogonal) first- 
order design previously discussed, is also rotatable. 

To obtain the conditions for rotatability in a second-order design, it better serves 
the purpose of this paper if most of the theory is omitted and emphasis is placed on only 
essential features. (For detailed information, see refs. 2 and 5 to 7.) In the present 
paper, it is necessary to define the following quantities in terms of the normalized vari- 
ables Xj, Xg, . . ., x k *. 



Notice that all of these quantities are elements of the matrix X’X. Moreover, p. is the 
fourth pure moment while if the fourth mixed moment. In terms of p^ and m^, 

the conditions for rotatability in a second-order response surface are as follows: 

(1) All odd moments through order four equal 0 


(2) The ratio of the fourth pure moment to the fourth mixed moment equals 3 (that 

P< \ v 


. k; i * j 



n 

n 

n 

1 V „ 

1 V 

1 v 

n Z x ih’ 

n Z X ih jh’ 

n Z 

h=l 

h=l 

h=l 


condition U; requires mat moments sucn as - ^ ^ x ih x jh ? n £j x ih x jh’ 

n n k=l h=l h=l 

n Z X ih ’ n Z x i h 3 x jh , and so forth, which are elements of the matrix X*X, all 


h=l h=l 

equal 0 for all i and j when i * j. 


Consider imposing rotatability on a central-composite design. By the nature of this 
design, it is easy to verify that condition (1) is satisfied. Thus, in a central-composite 
design, rotatability is equivalent to determining the value a which satisfies the require- 
ment that 
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By using the matrix X given in equation (14) as an example, it can be generalized that 


n 


1 x ih 


4 = 2 k + 2 a 4 


h=l 


and that 
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2 x 4 J = 2 k 


h=l 


ih A jh 


for any central -composite design. Therefore , the value of a for which 


2 k + 2a 4 


= 3 


is 


a = 2 


k/4 


(18) 


This value assures rotatability in a central-composite design for a given k. 

To illustrate a rotatable second-order central -composite design, consider the case 
of k = 2 controlled variables. By using equation (18), it is readily determined that 
a = for k = 2. Thus, the geometric representation of the desired region is the 
octagon of the following sketch: 



The design points are at the center and at the eight vertices of the octagon. Therefore, 
the design matrix X follows from this configuration and is given by 
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x 0 

X 1 

x 2 

x l 2 

x 2 2 

X 1 X 2 

_ 1 

-1 

-1 

1 

1 

1 

1 

-1 

1 

1 

1 

-1 

1 

1 

-1 

1 

1 

-1 

1 

1 

1 

1 

1 

1 

1 

0 

0 

0 

0 

0 

1 


0 

2 

0 

0 

1 

if 

0 

2 

0 

0 

1 

0 

-1/2 

0 

2 

0 

1 

0 

12 

0 

2 

0_ 


12 4 

It can be verified that all odd moments through order four equal 0, , m^ = -, and 

rotatability is thus attained. Perhaps it is important to remark here that to convert nor- 
malized units equal to ±\/2 or 0 to corresponding natural units, equation (2) can be used to 
solve for the natural variable. 

In addition to rotatability, another equally important property is now considered for 
a second-order central-composite design. Recall that, by the definition of rotatability, 
the variance of a predicted response at any arbitrary point interior to the experimenta- 
tion region is a function of the distance r from the point to the center of the region. 
Since the variance increases as the distance from the center increases, Box and Hunter 
(ref. 6) develop criteria so that the error of a predicted response at r = 1 is the same 
as the error at the design center. Such a consideration is based on the notion that the 
prediction of a response anywhere inside the region for which r < 1 is of uniform 
importance and, therefore, the error should be the same. This attractive property of 
uniform error can be achieved through a judicious choice of the value for the fourth 
mixed moment m^ without loss of rotatability. Specifically, the value of which 

gives rise to a uniform error for predicting a response anywhere inside a region bounded 
by r = 1 is made possible by measuring the response more than once at the center of 
the design. For rotatable second-order designs, these values of m^ are given for up 
to and including six controlled variables as follows: 


■ 

k 

m ij 

2 

0.3187 

3 

.4093 

4 

.5106 

5 

.6120 

6 

,7056 


15 




For a central-composite design, there remains to derive an expression which, for 
a given k, relates the uniform -error value of m^ to the corresponding number of mea- 
surements of the response needed at the center of the design. Let the total number of 
measurements be 


n = 2 k + 2k + n c 

where n c indicates the number of measurements at the center. Since m^ is the dot 
product of the squares of the ith and jth columns of X divided by n, then for a central- 
composite design 


2 k 

m ij = "k 

J 2 k + 2k + n c 

In solving equation (19) for n e , the relation reduces to 


(19) 


n c = ^-(2 k + 2k) (20) 

™ij 

where the value of m^ has already been defined for k § 6 controlled variables. Since 
n c must be an integer, values obtained by equation (20) have to be rounded off to the near- 
est whole number and, thus, absolutely uniform variance cannot be achieved. The neces- 
sary parameters for a second-order rotatable central-composite design with nearly uni- 
form variance are given, for up to six controlled variables, as 


k 

Value of a 
for 

rotatability 

No. of measurements 
at center 

for uniform error 

Total no. 
of measurements 

2 

1.414 

5 

13 

3 

1.682 

6 

20 

4 

2.000 

7 

31 

5 

2.378 

10 

52 

6 

2.828 

15 

91 


ANALYSIS OF SECOND-ORDER PREDICTION EQUATIONS 

Although the primary task of this report is the discussion of optimal -design con- 
cepts for the determination of a prediction equation, the task is not complete without at 
least a brief discussion of the analysis of such an equation, especially with respect to 
locating local maxima or minima. Once the prediction equation of a second-order 
response surface is available, it is usually desirable to determine the coordinates of the 
controlled variables at which the response is maximum (or minimum). In locating this 
point, let the predicted response y at any point s' = j~s^ Sg . 
terms of the normalized variables as 


s kJ be written in 
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( 21 ) 


where 


and 


y = 0 q + s’y + s'Bs 


>.t 

y = 

^2 * * ' 

^k) 


r— 

n 

^ / 


%i hi/ 2 

• • • hiP 


$n/ 2 %2 

• • • hi/ 2 

B = 

• • 

• • 

• • 
• • 


• • 
^lk/ 2 ’ * • 

. • 

A 

^kk 


( 22 ) 


Let Sg be the coordinates of the point at which the response is maximum (or minimum) 
within the experimentation region. The vector s Q is found by solving the equation 

d £=0 

ds 

Thus, by differentiating equation (21) and by virtue of the symmetry in B, 

££ = y + 2Bs 
ds - 


which, upon simplification, reduces to the needed relation 

5o = -K^ 


(23) 


The predicted response y Q corresponding to s Q is determined, by substituting equa- 
tion (23) into equation (21) and simplifying, as 

* 1 1 
y 0 = ^0 _ 2 rB ^ + 4 — B y - 




= Po + 2?bz 


(24) 


Therefore, if in fact there exists a point within the experimentation region which yields 
maximum (or minimum) response, the coordinates of such a point are given by equa- 
tion (23) while the maximum (or minimum) response is given by equation (24). The reader 
is cautioned that the existence of either a maximum or a minimum at Sq is not guaran- 
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teed because of the possibility that Sq may be a saddle point. Thus, it is imperative to 
determine the nature of the response at the point Sq. In pursuit of this goal, con- 
sider rewriting equation (21) by translating it from the current origin (design center) to 
the point Sq. This involves defining a new vector t_ = s - s^. As a result, equation (21) 
becomes 

y = j3 0 + L’r + sjj y + (t_’B + gJ,B^ + Sq) 

= ^0 + + io B ®0 + I'z + L’BSp + + t^Bt_ 

= 0 O + s ' Q y + £qBs q + t_’(y + 2 BSq) + t*B^ 

However, 0Q + s^y + SqBSq is the predicted response at the point Sq, and y + 2 BSq = 0. 

Therefore, the translation of equation (21) from the design center to the point Sq results 
in the relation 

y = y 0 ( 25 ) 

Since B, as defined by equation (22), is a real symmetric matrix, the principal-axis 
theorem is invoked to rotate equation (25) to a new set of variables w^, w 2 , . . ., w^ 
and thus reduce it to its canonical form, from which the nature of the response at _Sq is 
immediately apparent. The rotation is accomplished by an orthogonal transformation of 
the form 

t = Mw (26) 


where w* = jw^ w 2 . . . w^j and M is the orthogonal matrix of the transformation. 
Substituting equation (26) into equation (25) yields 

y = y 0 +w ,M, B M w (27) 


However, since B 


is symmetric, then 


m'bm 


0 0 
X 2 0 


0 

0 


where X^, X 2 , 


|_° 0 

., A^ are the eigenvalues of B. Hence, 
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and the relation given by equation (25) reduces to its canonical form given by 

k 

y = y 0 + Y, x i w i 2 ( 28 ) 

i=l 

Equation (28) was derived from equation (21) by translating the latter to Sq and then 
rotating it to the principal axes w^, W£, . , ., w^. As a result, the nature of the response 
at the point Sq is manifested by determining the eigenvalues Xj, Xg, * • •> of B. 

Through an inspection of the canonical form given by equation (28), the interpreta- 
tion of the response at Sq is as follows: 

(1) The response £q is maximum if all eigenvalues are negative 

(2) The response £q is minimum if all eigenvalues are positive 

(3) The response is neither maximum nor minimum if all of the eigenvalues 

are not of the same sign; this constitutes a saddle point at Sq 

If the response at Sq is maximum, a move in any direction away from Sq results in a 
decrease of the response. Conversely, if the response at Sq is minimum, a move away 
from Sq produces an increase in the response. In either case, the more pronounced 
change occurs along that axis whose corresponding eigenvalue has the greatest magnitude. 
When a saddle point is encountered, a move away from Sq along an axis w^ may result 
in an increase or a decrease of the response depending on whether the ith eigenvalue is 
positive or negative. Finally, it is necessary to remark that it is possible for the coor- 
dinates of the point Sq to fall outside the region under consideration. If such is the 
case, further experimentation in the direction of the point Sq is recommended depending, 
of course, on the experimenter's goals. 

Again, it is beneficial to provide an example to illustrate the procedures for a 
second-order model. Consider the situation in which a researcher attempts to determine 
criteria for regenerative adsorbents appropriate for CO 2 control in manned spacecraft. 

Of specific interest is the ability to predict CO 2 adsorption rate across a gas film layer 
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as a function of velocity V, CO 2 partial pressure P, and plate temperature T. To 
determine a second-order prediction model, 20 measurements are taken according to a 
rotatable central-composite design where the response is the total grams of CO 2 
adsorbed by the plate in a 10 -minute period. ' 

The relations between natural and normalized variables are 

x _ V - 45.72 

1 22.86 


Moreover, the range of experimental investigation for each of the three controlled vari- 
ables is given in both normalized and natural units as follows: 


Normalized 

units 

Velocity, 

m/min 

CO 2 partial 
pressure, 
mm Hg 

Plate 

temperature, 

°C 

-1.682 

7.32 

3.1 

0.54 

- 1.000 

22.86 

6.7 

6 

.000 

45.72 ! 

12.0 

12 

+ 1.000 

68.58 

17.3 

22 

+1.682 

84.13 

20.9 

27.46 


The X and x'x matrices, the vector of measurements y, and the vector X’y 
for this example are given as follows: 






2 

2 

2 




X G 

X 1 

X 2 

X 3 

X 1 

X 2 

X 3 

X 1 X 2 

X l x 3 

X 2 X 3 

1 

-1 

-1 

-1 

1 

1 

1 

1 

1 

1 

1 

1 

-1 

-1 

1 

1 

1 

-1 

-1 

1 

1 

-1 

1 

-1 

1 

1 

1 

-1 

1 

-1 

I 

t 

1 

-1 

1 

1 

1 

1 

-1 

-1 

1 

-1 

-1 

1 

1 

1 

1 

1 

-1 

-1 

1 

1 

-1 

1 

1 

1 

1 

-1 

1 

-1 

1 

-1 

1 

1 

1 

1 

1 

-1 

-1 

1 

1 

1 

1 

1 

1 

1 

1 

1 

1 

1 

1 

- 1.682 

0 

0 

2.828 

0 

0 

0 

0 

0 

1 

1.682 

0 

0 

2.828 

0 

0 

0 

0 

0 

1 

0 

- 1.682 

0 

0 

2.828 

0 

0 

0 

0 

1 

0 

1.682 

0 

0 

2.828 

0 

0 

0 

0 

1 

0 

0 

- 1.682 

0 

0 

2.828 

0 

0 

0 

1 

0 

0 

1.682 

0 

0 

2.828 

0 

0 

0 

1 

0 

0 

0 

0 

0 

0 

0 

0 

0 

1 

0 

0 

0 

0 

0 

0 

0 

0 

0 

1 

0 

0 , 

0 

0 

0 

0 

0 

0 

0 

1 

0 

0 

0 

0 

0 

0 

0 

0 

0 

1 

0 

0 

0 

0 

0 

0 

0 

0 

0 

1 

0 

0 

0 

0 

0 

0 

0 

0 

0 
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x’x = 


20 

0 

0 

0 

13.656 

13.656 

13.656 

0 

0 

0 


0 

13.656 

0 

0 

0 

0 

0 

0 

0 

0 


0 0 

0 0 

13.656 0 

0 13.656 

0 0 

0 0 

0 0 

0 0 

0 0 

0 0 


13.656 13.656 

0 0 

0 0 

0 0 

.24 8 

8 24 

8 8 

0 0 

0 0 

0 0 


13.656 0 0 0 

0 0 0 0 

0 0 0 0 

0 0 0 0 

8 0 0 0 

8 0 0 0 

24 0 0 0 

0 8 0 0 

0 0 8 0 

0 0 0 8 


1.996 

2.812 

5.171 

5.035 

2.676 

4.173 

3.266 

6.668 

2.359 

3.583 
2.449 
5.851 
2.223 
2.994 
6.519 
6.218 
6.354 

5.583 
6.491 

J 5.781 


x r y- 


88.2020 

7.6377 

14.2052 

3.0658 

48.6010 

55.2694 

46.5507 

0.9530 

4.2190 

- 2.3130 


Solving for the coefficients of the second-order model gives the prediction equation in 
terms of the normalized variables as 

y = 6.1283 + 0.5596X! + 1.0406x 2 + 0.2246x 3 - 0.9350XJ 2 - 0.5181x 2 2 

- 1.0632x 3 2 + 0.1191x^2 + 0.5274x^2 - 0.2891x 2 Xg 

The prediction equation in terms of the natural variables reduces to 

y = -12.2808 + 0.13591V + 0.68956P + 0.44319T - 0.001789V 2 - 0.01844P 2 

- 0.016613T 2 + 0.000983VP + 0.002884VT - 0.006818PT 
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By appealing to an analysis of the second-order prediction equation, the matrix B 
as defined by equation (22) is found to be 


B = 


-0.9350 

0.0595 

0.2637 


0.0595 

-0.5181 

-0.1445 


0.2637 

-0.1445 

-1.0632 


while the vector is 

y* = jj).5596 1.0406 0.2246] 

Therefore, by using equation (23), the normalized coordinates of Sq are 


0.382 


£o ~ 


1.031 

0.060 


x 

X 

X 


1 

2 

3 


The coordinates are clearly inside the region of investigation. Moreover, the amount of 
CC>2 adsorbed at Sq is estimated to be Yq = 6.7780 grams. As the analysis is con- 
tinued, the eigenvalues of B are found to be Xj = -0.4820, X£ = -0.7347, and 
Xg = -1.3000 and thus the canonical form of the prediction equation is 

y = 6.7780 - 0.4820w 1 2 - Q.7347w 2 2 - 1.300Wg 2 (29) 

Since all eigenvalues are negative, the amount of CO2 adsorbed at s^ is maximum. In 
terms of natural units, therefore, maximum CO2 adsorption occurs at a velocity of 
54.45 m/min, a CO2 partial pressure of 17.5 mm Hg, and a plate temperature of 14.48° C. 
From equation (29), it is apparent that the amount of CO2 adsorbed, away from Sq, is 
more sensitive to change in either w 2 or Wg than in w^. 


CONCLUDING REMARKS 


Before scientific and engineering measurements are made for determining empirical 
relationships between variables of interest, it is always advantageous to devote consider- 
able effort to identifying such things as what is to be measured, where it should be mea- 
sured, and how much is to be measured. The validity of the derived empirical model 
ultimately depends on such considerations. Therefore, by considering two low -order 
polynomial models, the essence of this report is to bring forth criteria necessary for well- 
planned experimental tests. In fact, the techniques described in this report are not only 
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essentially optimal in the sense of yielding prediction equations which minimize the error 
between a measured and a predicted response but are also economical in that they require 
a relatively small number of measurements. 

Langley Research Center, 

National Aeronautics and Space Administration, 

Hampton, Va., October 26, 1972. 
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APPENDIX 


PROOF OF THEOREM OF EQUATION (8) 

The theorem of equation (8), which is to be proved, is stated as follows: If the 
design matrix X is selected such that each column of X is mutually orthogonal to 
every other column in X, then the first-order approximation to the response surface is 
optimal in the sense that the residual variance of a predicted response is minimized. 

Consider the determinant of X'X. Let a^ be the (h,h) element of X T X and let 
Ahh be the cofactor of a^ in X’X. Moreoever, let A^.^ be the cofactor of 
in A^ for i,j + h. Since X is a matrix of full rank m < n, the matrix X’X is a 
positive definite symmetric matrix; thus, consider Cauchy's expansion (ref. 8) of |x’x| 
given by 

m m 

l x ’ x l = a hh A hh " J ^ a hi a hj A ij:hh 

i=0 j=0 


- a hh A hh “ Q 


(h = 0, 1, 2, . . ., m; i,j * h) 


where a^A^ is a scalar and Q is a positive definite quadratic form in a^ for 
h * i. Multiplying through by a 2 and dividing by |x’x| yields 

a 2 = g2 ( a hh A hh) _ g 2 Q 

|x’x| jx'x| 

Since a^ is the hth diagonal element of x'x, then a^ is the sum of squares of the 
elements in the hth column of X. However, these elements are all either -1 or +1; thus, 


(Al) 


"hh 

a^v, = n. Moreover, - 

hh | x’x) 


is the hth diagonal element of (X'X) - *. Then, by definition, 


var 


a 2 A 


hh 


x’x! 


A A 

where (3^ is the hth component of the estimate vector £, Substituting into equation (Al) 
yields 


cr = n var 




or 


m(U = ¥ * - 4 f 1 + r 2 -) 

' h ' n| X x| n \ |x’x|j 
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APPENDIX - Concluded 


Since X r X is positive definite, Ix’xl > 0; moreover, since Q is positive definite, the 
term . ^ , can be no smaller than 0. In fact, it can equal 0 only when all the off- 

Ix’xl 

diagonal elements a hi equal 0 for h * i t This implies that x'x is a diagonal matrix 
which, in turn, implies that each column of X is mutually orthogonal with every other 
column in X. Therefore, var(^) is minimized when x’x is of the form given by 
equation (8). 
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